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ABSTRACT 

A low Reynolds number multiple-time-scale turbulence model (LMS) 
and its application to fully developed turbulent channel flows and 
pulsating pipe flows are presented. The LMS can describe the inequilibrium 
turbulence phenomena down to the viscous sublayer. The calculated fluid 
flow and turbulence fields for the channel flows are in better agreement 
with the direct numerical simulation (DNS) results than those obtained 
using a Reynolds stress turbulence model, and the calculated near-wall 
dissipation rates are in qualitatively correct agreement with the DNS 
results. The LMS also successfully predicts the rapidly varying phase-lead 
of the wall shearing stress that occurs in a narrow range of the 
dimensionless frequency (co + = cov/ux^) for the pulsating pipe flows while 

various other turbulence models fail to predict this phenomenon, and the 
LMS yields significantly improved numerical results for a wide range of the 
dimensionless frequency compared with those obtained using a rapid 
distortion theory (RDT). 


*NASA Resident Research Associate at Lewis Research Center. 



NOMENCLATURE 

Cpi model constants for e p equation (i=l,3) 

c t j model constants for e t equation (i=l,3) 

cp eddy viscosity coefficient 
Cpf constant coefficient (=0.09) 

f frequency 

fp wall damping function for eddy viscosity equation 

k turbulent kinetic energy (k=k p +k t ) 

k + normalized turbulent kinetic energy (=k/u x 2 ) 

k p turbulent kinetic energy in production range 

k p + normalized k p (=k p /u t 2 ) 

k{ turbulent kinetic energy in dissipation range 

k t + normalized k t (=k t /u T 2 ) 

p pressure 

P r production rate 

Re Reynolds number 

Ry turbulent Reynolds number (=fk y/v) 

T period of oscillation 

u b bulk velocity 

u time averaged velocity 

u + non-dimensional velocity (=u/u T ) 

u T friction velocity (=Vt w /p) 

u'v' Reynolds stress 

u'v ,+ normalized Reynolds stress (=uV/u x 2 ) 
xj spatial coordinates (={x,y,z}) 

y + wall coordinate (=yu x /v) 

8 S Stokes layer thickness (=V2 v/cd) 

e p energy transfer rate 
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£p + nomalized energy transfer rate (=vep/u T 4 ) 

e t dissipation rate 

e t + nomalized dissipation rate (=ve t /u T 4 ) 

p molecular viscosity 

v kinematic viscosity 

t w wall shearing stress 

Pt turbulent viscosity 

p density 

Gj turbulent Prandtl number for i-equation, i={kp,Ep,kt,£t} 

to angular velocity (=2rcf ) 

to + dimensionless frequency (=cov/u t 2 ) 

Superscripts 

fluctuating velocity component 
( ) time-averaged value 

( ) ensemble-averaged value 

l( )! amplitude of oscillation 

Subscripts 

j denotes spatial coordinates (={x,y,z}) 
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1. Introduction 

Either the development of a new near-wall turbulence model or the 
improvement of various near-wall turbulence models are hindered due to 
the lack of detailed knowledge on the near-wall turbulence structure 
especially that for the near-wall dissipation rate. Recently, Kim et al. 
(1987) provided detailed near-wall turbulence data for fully developed 
channel flows obtained from direct numerical simulations. The validity of 
the DNS results is established by comparing them with measurable 
turbulence quantity. A near-wall dissipation rate obtained from the DNS 
results and that proposed by Patel et al. (1985) are shown in figure 1. The 
dissipation rate obtained from the DNS results show that the peak value is 
located at the wall while the dissipation rate model proposed by Patel et al. 
(1985) shows that the dissipation rate attains its peak value at y + ~12. 
Numerical calculations using various turbulence models yield near-wall 
dissipation rates similar to the one proposed by Patel et al. (1985). It is not 
clear if there exist any turbulence model that can yield a near-wall 
dissipation rate that is in qualitatively correct agreement with the DNS 
result and that can also yield accurate numerical results for a few different 
classes of turbulent flows. Consequently, the influence that the two 
distinctly different near-wall dissipation rates can cause on numerical 
calculations of turbulent flows is not known yet. The discrepancy that can 
be caused by the different near-wall dissipation rates are discussed in this 
paper. 

Numerous efforts have been made to develop or to improve turbulence 
models to correctly predict steady simple shear layers as well as complex 
turbulent flows. However, only a very small amount of efforts has been 
made for numerical investigations of unsteady turbulent flows. Calculations 
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of a pulsating pipe flow using an algebraic turbulence model, a k-e 

turbulence model, and a Reynolds stress turbulence model (RSM) can be 
found in Kebede et al. (1985). It can be seen in the reference that the 
algebraic turbulence model yields the most accurate numerical results and 

the RSM yields the worst numerical results, and the cause for the 
deteriorated numerical results obtained using the RSM is not clearly 
known. Mao & Hanrathy (1986) carried out experimental and numerical 

investigations of pulsating pipe flows. Calculations of the pulsating pipe 

flows at a wide range of the dimensionless frequency (©+) revealed that an 
algebraic turbulence model discussed in the Thorseness et al. (Mao & 
Hanrathy, 1986) yields qualitatively inaccurate numerical results while a 
pressure-relaxation algebraic turbulence model yields numerical results 
that are in good agreement with the measured data. However, it is 
admitted in Mao & Hanrathy (1986) that the pressure-relaxation model 
lacks a theoretical background. More recently, Mankbadi & Liu (1992) 
developed an algebraic turbulence model based on the rapid distortion 
theory (RDT) and showed that the turbulence model yields accurate 
numerical results for pulsating pipe flows at high dimensionless frequency 
while the turbulence model still can not yield accurate numerical results 
for the flows at low dimensionless frequency. As is obvious from the above 
discussion, numerical investigations of unsteady turbulent flows are mostly 
made using algebraic turbulence models since k-e or RSM does not yield 
accurate numerical results. However, applicability of algebraic turbulence 
models is limited since these turbulence models can not yield accurate 
numerical results for various complex turbulent flows. 

The multiple-time-scale turbulence model presented in this paper is 
based on a simplified split-spectrum method (Hanjelic et al. 1980; Kim & 
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Chen 1989). In the high Reynolds number multiple-time-scale turbulence 
model presented in Kim & Chen (1989), the near-wall turbulence field is 
described using a wall function method. This turbulence model yields 
accurate numerical results for complex turbulent flows in case the 
turbulence field is not strongly influence by the near-wall turbulence 
structure. Such turbulent flow cases can be found in a turbulent flow over 
a backward -facing step, a wall-jet flow, and confined coaxial jets with and 
without swirling velocity components. Later, the predictive capability of 
the turbulence model was improved by incorporating a "partially low 
Reynolds number near-wall turbulence model" into the multiple-time-scale 
turbulence equations (PLMS). In the PLMS (Kim 1990, 1991; Kim & Benson 
1992), the near- wall turbulence field is described by a single time-scale 
turbulence model derived from a k-equation turbulence model. The PLMS 
yields significantly improved numerical results for the turbulence 
structure in the near-wall region. It can be found in Kim (1990) that the 
PLMS successfully predicts the extensive growth of the shock- separated 
recirculation region in the transonic flow over a curved hill. It can also be 
found in Kim & Benson (1992) that the PLMS can accurately resolve the 
strong interaction between circular jets and crossflows and successfully 
predicts the horse-shoe vortex located along the circumferences of the jet 
exits. Unfortunately, the turbulence model still can not yield accurate 
numerical results for a class of turbulent flows in which the entire fluid 
flow is strongly governed by the near-wall turbulence structure since, in 
such cases, a large portion of the turbulence field is described by a single 
time-scale k-equation turbulence model. This class of turbulent flows can 
be found in fully developed channel flows at very low Reynolds numbers 
and unsteady pulsating pipe flows. In this paper, a new improvement is 
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incorporated into the multiple-time-scale turbulence model to accurately 
resolve this class of turbulent flows and to extend the applicability of the 
turbulence model to wider classes of turbulent flows.. 

The time averaged or the ensemble averaged Navier-Stokes equations 
and the turbulence equations are solved by a finite volume method that 
incorporates a pressure-staggered mesh and a partial differential equation 
for incremental pressure (Kim 1990, 1991; Kim & Benson 1991, 1992). 
Various laminar and turbulent flows solved using the numerical method 
include: reattaching shear layers in a divergent channel, transonic 

turbulent flows with shock wave - boundary layer interaction, three- 
dimensional turbulent flows of circular jets in crossflows, and self- 
sustained unsteady flows over a circular cylinder and a square cylinder. 
Details of the numerical method and the independence tests of the 
numerical method on the grid size and on the time-step size can be found 
in Kim & Benson (1991) and the references cited therein. It can also be 
found in these references that the present numerical method yields 
accurate numerical results for separated and recirculating flows using 
highly graded meshes. 

2. Low Re multiple-time-scale turbulence equations 

In the multiple -time-scale turbulence models, the turbulent transport of 
mass and momentum is described using the time-scale of large eddies and 
the dissipation rate is described using the time-scale of fine-scale eddies 
(Kim & Chen 1989; Kim & Benson 1992). In this regard, the convection- 
diffusion equations of the multiple-time-scale turbulence models describe 
the physically observed turbulence phenomena most naturally. The use of 
multiple time-scales also enables the turbulence equations to resolve the 
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inequilibrium turbulence" phenomena and to model the cascade process of 
the turbulent kinetic energy. Here, the "inequilibrium turbulence" 
represents the state of a turbulence field in which P r /e t varies rapidly in 

space so that the shape and the frequency domain of the spectral density 
varies widely in space. The influence of the inequilibrium turbulence on 
turbulent transport of mass and momentum can be observed in a number 
of semi-emperical data (theoretically derived data from a set of measured 
data). Detailed analyses of these capabilities can be found in Kim & Benson 
(1992) and the pertinent results are presented in the "Inherent capability 
of LMS equations" sub-section. 

It can be seen in Kim & Benson (1992) and the references cited therein 
that the numerical results for various complex turbulent flows (e.g., 
turbulent flows subjected to extra strains caused by streamline curvature, 
interaction of multiple number of shear layers, and shock wave-boundary 
layer interactions) obtained using the multiple-time-scale turbulence 
model are in as good agreement with the measured data as those obtained 
using an optimized k-e turbulence model, algebraic Reynolds stress 
turbulence model (ARSM) or RSM for each flow. Undoubtedly, the 
anisotropy of the turbulence is the most easily detectable phenomenon in a 
measurement of a turbulent flow and hence, in theoretical investigations of 
turbulence, the emphasis was laid upon improving the ARSM and the RSM. 
However, numerous numerical investigations carried out during the last 
decades show that the ARSM and RSM still can not accurately predict the 
turbulence phenomena occurring in various flows unless the pressure- 
strain rate correlation is optimized for each flow. The capability of the 
multiple-time-scale turbulence model to solve widely different complex 
turbulent flows is attributed to its capability to resolve the inequilibrium 
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turbulence and to model the cascade of the turbulent kinetic energy. A low 
Reynolds number multiple-time-scale turbulence model that can resolve 
the inequilibrium turbulence phenomena down to the viscous sublayer and 
calculations of turbulent flows that are highly sensitive to the near-wall 
turbulence structure are presented in this paper. 


2.1. LMS equations 

The multiple-time-scale turbulence equations that include the near-wall 
modifications are described below. The turbulent kinetic energy (kp) and 
the energy transfer rate (Ep) equations for energy containing large eddies 

are given as; 


^ P kp) + A ( p Ujkp) 





= pP r - P e p 


( 1 ) 


l(pe p ) + A(pu jEp ) 



where p is the density, uj (={u,v}) is the ensemble-averaged velocity, p is 
the molecular viscosity, p^ * s the turbulent viscosity, oiq) and o £ p are the 
turbulent Prandtl numbers for the kp- and Ep -equation, respectively, Cpj 
(i=l,3) are model constants, and the production rate (P r ) is given as 


P.=^ 2 



3u + 3vj^\ 

idy dJ I 


The turbulent kinetic energy (k t ) and the dissipation rate (e t ) equations for 
fine scale eddies are given as: 


d , , . 3 . , . a f p t 3k t 

— (pk t ) + — (pujk t ) - — p + 1 

3t 3xj 3xj |v c kt' 3x| 


= pe p - pe t 


(3) 
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where o k t and C7 £ t are the turbulent Prandtl numbers for the kt* and ep 
equation, respectively, and c^j (i=l,3) are model constants. The turbulent 
eddy viscosity is given as; 

f k 2 

^ = P c pf f u-^ (5) 

The turbulence model constants for the LMS equations are given as; 
c pi (i=1 ’ 3)= (0*21, 1.32, 1.84}, c ti (i=l,3)= {0.32, 1.21, 1.65}, a kp = a kt = 0.75, 

°ep = a et = 1*15, and c p f = 0.09. A few model constants in c p j and c t j are 

slightly different from those used in the PLMS. Calculations of free shear 
layers and complex turbulent flows using the LMS show that the 
inequilibrium turbulence structure in the external region is not influenced 
significantly by the slightly re-distributed model constants. Furthermore, 
many complex turbulent flows are not strongly influenced by the near-wall 
turbulence structure so that the numerical results for complex turbulent 
flows obtained using the LMS are only slightly better than or almost the 
same as those obtained using the PLMS. 

The high Reynolds number e p - and epeauation can not correctly 
describe the inequilibrium turbulence phenomena in the near-wall region 
since the wall proximity and the molecular viscosity dominate the 

development of the turbulence field in the region. The wall damping 
functions f £p and f £t are constructed in such a way that the functions 

reproduce physically consistent near-wall distributions of e p and e t and 

render the partial differential equations well posed at y~0. The wall 


damping functions f n> f ep and fet are described later in the "Near-wall 
turbulence" sub-section. 


2.2. Near-wall e p and e t 

The near-wall £p and e^ are obtained from analytical solutions of the kp- 
and k t - equation. Adding up equations (1) and (3) yields 


A ( p k) + g. (pu . k) .l.j(, + a)^j = p p r . p e 1 


S I 


( 6 ) 


where k (=kp+k[) is the turbulent kinetic energy. At very close to the wall 

(y~0), the molecular diffusion term and the dissipation rate dominate over 
the other terms in equation (6) since k=0, Pr=0 and p l =0. Formally 
integrating -3 { jo(3 k/5 x j )}/3x j = -pe t yields 

e t = 2vk/y 2 (7) 


where v=p/p is the kinematic viscosity. The dissipation rate obtained using 
equation (7) and the turbulent kinetic energy of DNS results for Re(u x ) = 
180 are compared with that of DNS result in figure 1, where u x (=Vx w /p) is 
the friction velocity and x w is the wall shearing stress. The near-wall 
dissipation rate obtained using e t = 2v(dfk"/0y) 2 is also shown in the figure. 
The latter expression is used in a number of low Reynolds number 
turbulence models cited in Patel et al. (1985) and in the RSM proposed by 
Lai & So (1990). It can be seen in the figure that the near-wall dissipation 
rate obtained using equation (7) compares more favorably with the DNS 
result for a wider range of y + than that obtained using e t = 2v(9lfk"/9y) , 
where y + =yu x /v is the normalized distance from the wall. It can also be 

seen in the figure that the slope of the dissipation rate obtained using 
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equation (7) is less steeper than that of the DNS result. The less steeper 
slope is caused by neglecting the turbulent diffusion term in deriving 
equation (7). Similarly, the energy transfer rate for y~0 is obtained by 
formally integrating equation (2) and is given as 

e p = 2vk p /y 2 (8) 

In the present calculations, the near-wall e p and e t for R y <5 (which 
corresponds to y + <2~3) are obtained using equations (7) and (8), where 
^y = ^ k "y/ v is the turbulent Reynolds number based on the distance from 
the wall. 


2.3. Near-wall turbulence 

The instantaneous velocities in the viscous sublayer can be written as 
(Patel et al. 1985) 


u' = ajy + a 2 y 2 + 

v’ = b 2 y 2 + 

w’ = qy + c 2 y 2 + 


} 


(9) 


where aj, bj, and cj are functions of time, aj = b- = cj = 0, and the over-bar 

denotes a time-averaged value. The fluctuating normal velocity grows in 
proportion to the square of the distance from the wall due to the wall 
proximity. Thus the turbulent kinetic energy and the Reynolds stress in the 
region very close to the wall (y~0) can be written as; 

k = ^- (u' 2 + v' 2 + w' 2 ) = -L (a 2 + c 2 ) y 2 

— 2 — , 2 1 ( 10 ) 
-u'v' = - a,b 2 y 3 J 
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Equation (10) indicate that the growth rates of the turbulent kinetic energy 
and the Reynolds stress are proportional to the second and the third power 
of the distance from the wall, respectively. 

In the Boussinesq eddy viscosity assumption, the Reynolds stress is 
given as; 


-uV = v t 


8u 

By 


( 11 ) 


Solving equations (5) and (11) for f^ yields 


f .= 


-u v 


CV (k z /e p )Ou/9y) 


( 12 ) 


Consider e p = 2vk p /y = constant and u + =y + for y~0, where u + =u/ux- Hence, 


f M ~ - — > oo for y — > 0 

u y j 


(13) 


Inside the viscous sublayer (0<y + <5), the molecular viscosity dominates 
over the turbulent viscosity; and in the fully turbulent region (y + >30~40), 
the turbulent viscosity dominates over the molecular viscosity. Therefore, 
the wall damping function needs to satisfy the following constraints to 
correctly describe the physically observed near-wall turbulence field, that 
is 


f^«lfory + = 5 1 

f p -»l for y + 100 J 

The wall damping function in the LMS is given as 


( 14 ) 
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1 - exp(- Pi Vr^ - (3 2 R y - P 3 Ry) 
1 - exp(-p 4 R y ) 


( 15 ) 




and is shown in figure 2. The Taylor series expansion of equation (15) for 
y~0 clearly shows that f p - Pj/^Vr^) « 1/y for R y ~0. The coefficients 

Pi(i=l,4) = {0.005, 0.001, 0.00011, 0.14} are obtained by optimizing the 

constants to yield the best numerical results for the fully developed 
channel flow at Re(u x ) = 180. 

For y~0, all the terms except the one containing the ep in equation (2) 

vanish. For the partial differential equation to be well posed at y~0, f ep 

should vanish at y~0. Thus the wall damping function for the Ep -equation is 
given as; 

fep = 1 - a ep exp(-R y ) (16) 

where a e p = 1.0. For the Ej-equation, equation (4), all the terms except the 

load functions become negligible at y~0. Thus equation (4) becomes an 
algebraic equation governing the ratio of Et/ep at y~0. Since the production 
rate vanishes on the wall, the ratio of Et/£p needs to be slightly greater 
than unity at the wall for the dissipation rate equation to be consistent 
with the inequilibrium analysis described in the following sub-section. 
Thus the wall damping function is given as 

ht = 1 - a et exp(-R y ) (17) 

where a e t = 0.13. The use of Ry in the wall damping functions (fp, f ep and 
f et ) makes the influence of the wall proximity to disappear as the fully 

turbulent external region is approached and hence it renders the 
turbulence equations to be consistent with a physical observation that the 
turbulence length scale in the near-wall region is related to the normal 
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distance from the wall while that of the external flows is related to the 
flow field characteristics (Roshko 1976). 

2.4. Inherent capability of LMS equations 
The spectral density curves that correspond to different levels of P r /e t 

and a partition representing the simplified split-spectrum are shown in 
figure 3(a). The spectral density curves are constructed based on the 
measured data of Klebanoff (1955) and Wygnanski & Fiedler (1969). The 
spectral density curves show that the energy-containing large eddies 
generated by the instability of the mean fluid motion are characterized by 
low frequency and the fine scale eddies in the dissipation range are 
characterized by high frequency. The curves also show that the ratio of 
kp/kt is determined by the shape and the frequency domain of each 
spectral density curve and that the ratio of k p /k t is large for large eddies 
and is small for fine scale eddies. Thus the cascade of turbulent kinetic 
energy is characterized by the ratio of k p /k t . The capability of the 

multiple-time-scale turbulence equations to model the cascade process is 
achieved by solving the splitted turbulent kinetic energy equations (k p - 
and k t -equation). 

The inequilibrium analysis is based on the turbulence statics observed 
in the uniformly sheared flows (Tavoularis & Karnik 1989) and on the 
semi-empirical data used in the generalized algebraic stress turbulence 
models (Rodi 1972; Launder 1982; Kim & Chen 1988). For uniformly 
sheared flows and free-stream flows, the diffusion terms vanish. In such 
cases, dividing equation (1) by equation (3) yields (Kim & Benson 1992) 



The existence of asymptotic states in which P r /e t takes constant values that 

depend on the mean flow strain rates can be found in Tavoularis & Karnik 
(1989), and the dependence of e p /e t on the ratio of P r /e t can be observed 

in the semi-empirical data used in the generalized algebraic stress 
turbulence models as shown in figure 3(b). The eddy viscosity equation, 
equation (5), can be written in a form comparable with those of the 
generalized algebraic stress turbulence models and is given as 

. k 2 

= PS f n T- (19) 

where c^c^f (e t /e p ). Thus, in the multiple-time-scale turbulence model, the 
dependence of the eddy viscosity coefficient (c p ) on the inequilibrium 
turbulence (P r /et) is described by the ratio of e t /e p . The turbulence model 
constants (c p j and c t j) obtained by solving an under-determined system of 

equations that governs the growth rate of the turbulence intensity in the 

uniformly sheared flows and the decay rate of the grid turbulence produce 
a c p -curve shown in figure 3(b). The under-determined model constants, 

subjected to constraint conditions obtained from a near-wall equilibrium 
analysis (P r =e p = Et), are numerically optimized to yield best numerical 

results for simple shear layers such as a plane jet in a uniform stream and 
fully developed channel and pipe flows (Kim & Chen 1989; Kim 1991). It 
can be seen in figure 3(b) that e t /e p <l and the ratio of P r /e t increases 
faster than e p /e t (i.e., the slope of |(e t /e p )/(P r /e t )| is less than unity) for the 
production dominated case (P r /E t >l). In this case, P r >e p >e t so that equation 
(18) always yields a positive ratio of k p /k t , and the ratio of k p /k t is further 
increased as Pf/et * s increased. For turbulent flows in equilibrium state 
( p r =e p =e t)> equation (18) becomes indeterminate, and the ratio of k p /k t for 


16 



P r~£p®£t is determined by the imbedded constraint conditions. As the 
production rate becomes negligible (P r =0), P r <e t and the ratio of e p /e t 
decreases faster than P r /e t . In the latter case, P r <Cp<e t so that equation (18) 
always yields a positive ratio of k p /k t , and the ratio of k p /k t is further 
decreased as P r /et is decreased. Therefore the trend of k p /k t distribution 

obtained in the inequilibrium analysis is in correct agreement with that 
observed in the spectral density curves shown in figure 3(a). 

3. Computational results 

3.1. Fully developed channel flows 

Fully developed channel flows at Re(u x ) = 180 and 395 are considered 

below. In each case, the flow domain in x-coordinate direction extends 70 
channel heights in the downstream direction and that in the y-coordinate 
direction extends from the wall to the symmetric half of the flow domain. 
The computational domain is discretized by 97x91 grid points in x- and in 
y-coordinate directions, respectively. The grid size in the x-ccordinate 
direction is almost uniform while that in the y-coordinate direction is 
stretched by a factor of approximately 1.1. The smallest mesh size in the 
near- wall region is Ay + =0.25. The inlet boundary condition is obtained by 
appropriately scaling measured data for a fully developed channel flow so 
that the calculated flow field may reach a fully developed state at a shorter 
downstream location from the inlet boundary. The numerical results 
presented herein are obtained at 5.5 channel heights upstream of the exit 
boundary so that the numerical results are not obscured by the exit 
boundary condition either. Examination of numerical results obtained using 
a few different mesh densities and extents of the flow domain shows that 
the numerical results presented herein are grid independent. 
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The numerical results for the normalized mean velocity (u + =u/u x ), the 
turbulent kinetic energy (k + =k/u^), the Reynolds stress (-u'v' + = - u'v'/u 2 ) 
and the dissipation rate (e t +=ve t /u£) for Re(u t )=180 are compared with the 
DNS results and the numerical results obtained using a RSM (Lai & So 
1990) in figure 4. It can be seen in figure 4(a) that the calculated mean 
velocity profile is in excellent agreement with that obtained from the DNS. 
The growth rate and the peak value of the turbulent kinetic energy in the 
near-wall region obtained using the LMS is in better agreement with the 
DNS result than those obtained using the RSM as shown in figure 4(b). The 
calculated Reynolds stresses in figure 4(c) also show that the LMS yields 
improved numerical result than the RSM. The calculated dissipation rate is 
shown in figure 4(d). The slope of the dissipation rate obtained using the 
LMS is less steeper than that obtained from the DNS. The less steeper slope 
is caused by neglecting the turbulent diffusion term in deriving equation 
(7) and by the nonlinearity of the turbulence equations. The most 
distinguished difference between the LMS and the RSM as well as many 
other turbulence models can be observed in the calculated near-wall 
dissipation rates as shown in the figure. Extending the capability of the 
multiple-time-scale turbulence model to resolve the inequilibrium 
turbulence inside the near-wall layer enables the turbulence equations to 
yield a near-wall dissipation rate that is in qualitatively correct agreement 
with the DNS result. 

Detailed comparisons of the various terms in the budget of the turbulent 
kinetic energy are shown in figure 5, where e p + = ve p /u^ is the normalized 
energy transfer rate. It can be seen in the figure that each component 
obtained using the LMS overlaps with that obtained from the DNS in most 
part of the flow domain. The excellent comparisons between the LMS and 
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the DNS results indicate that the LMS can accurately resolve details of the 
near- wall turbulence structure. It can also be found in Lai & So (1990) that 
the present numerical results are in much better agreement with the DNS 
results than those obtained using the RSM, and the excellent agreement is 
attributed to the calculated near-wall dissipation rate that is in 
qualitatively correct agreement with the DNS result. The terms in the 
budget of kp- and k t -equation are shown in figures 5(b) and (c), 

respectively. It can be seen in the figures that the dissipation rate in the 
region very close to the wall is mostly balanced by the laminar viscous 
diffusion of the turbulent kinetic energy. The laminar viscous diffusion of 
the turbulent kinetic energy do not require any modelling and, hence, 

obtaining the correct dissipation rate at the wall depends on obtaining the 
correct near-wall distribution of the turbulent kinetic energy. Thus the 

qualitatively different near-wall dissipation rate obtained using the RSM is 
caused by the near-wall turbulent kinetic energy distribution that do not 
compare very well with the DNS result as shown in figure 4(b). Note that 
the near-wall energy transfer rate (ep) shown in figure 5(b) exhibits 

similar distribution as the dissipation rate obtained using single-time-scale 

% 

turbulence models in the sense that the peak value is not located at the 
wall. Thus the dissipation rate that attains its peak value at the wall is 
caused by a small amount of the laminar viscous diffusion of the 

dissipation range turbulent kinetic (k t ) at the wall as shown in figure 5(c). 

The composition of the turbulent kinetic energy is shown in figure 6(a). 
Since the peak production rate occurs at y + =12 and the balance of the 
turbulent kinetic energy budget in the region below y + = 12 is mostly 
achieved by the laminar viscous diffusion, the turbulent kinetic energy in 
the region is mostly contained in the production range. The ratio of kp/kj is 
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also shown in figure 6(a). The viscous sublayer region is entirely 
dominated by the viscous diffusion and, hence, the inequilibrium analysis 
developed based on the turbulence statistics of uniformly sheared flows do 
not apply exactly in the region. Away from the viscous sublayer, the ratio 
of k p /k t is decreased as P r /e t is decreased and, thus, the distribution of 
kp/k t is in correct agreement with the inequilibrium analysis. The near- 
wall distribution of e t /e p is shown in figure 6(b). It can be seen in the 
figure that e t /ep>l at the wall where the production rate vanishes, e t /£p< 1 
near y+=12 where the production rate attains its peak value and P r /e t >l, 
and e t /e p >l again in the region remote from the wall where P r /e t <l. Thus 
the calculated e t /e p is in correct agreement with the inequilibrium analysis. 

The numerical results for the mean velocity, the turbulent kinetic 
energy, the Reynolds stress and the dissipation rate for Re(u t )=395 are 

compared with the DNS results in figure 7. It can be seen in the figure that 
the calculated mean velocity profile, the turbulent kinetic energy and the 
Reynolds stress are in excellent agreement with those obtained from the 
DNS. Each term in the budget of the turbulent kinetic energy for 
Re(u x )=395 case exhibits improved comparison with the DNS results than 
the Re(u x )=180 case, and the improved comparison is caused by the 

turbulent kinetic energy and the near-wall dissipation rate profiles that 
are in closer agreement with the DNS results than the Re(u x )=I80 case. 

Measured data show that the growth rate and the peak value of the 
turbulent kinetic energy are increased as the Reynolds number is increased 
(Patel et al. 1985). The turbulent kinetic energy obtained from the DNS also 

show that the growth rate and the peak value of the turbulent kinetic 
energy for Re(u x )=395 are larger than those for Re(u x )=180. Since the near- 

wall dissipation rate depends mostly on the growth rate and the peak 
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value of the turbulent kinetic energy near the wall, the peak value of the 
dissipation rate at the wall for Re(u t )=395 is substantially larger than that 

for Re(u x )=180. It can be seen in figure 7(d) that the LMS correctly predicts 
the increased near-wall dissipation rate. 

3.2. Pulsating pipe flows 

Pulsating pipe flows at a wide range of dimensionless frequency 
(g> + = 0.0075~0.21) are considered below. In each case, the flow domain in x- 
coordinate direction extends 70 pipe diameters in the downstream 
direction and that in the y-coordinate direction extends from the wall to 
the center of the pipe. The computational domain is discretized by 144x64 
grid points in x- and in y-coordinate directions, respectively. The grid size 
in the x-ccordinate direction is almost uniform while that in the y- 
coordinate direction is stretched by a factor of approximately 1.1. The 
smallest mesh size in the near-wall region is approximately equal to 1/15 
5 S (which corresponds to Ay + =0.25~0.8 for different Reynolds numbers), 
where 5s=V2v/co is the Stokes layer thickness, a>=27tf is the angular velocity, 
and f is the frequency. This mesh is fine enough to resolve the influence of 
the imposed frequency on the near-wall fluid flow and to yield a grid 
independent solution. The transient bulk velocity is prescribed as 
u b = u^( l-|u^j sin(cot)}, where () denotes the time-averaged value, () 
denotes the ensemble-averaged value, and |( )l denotes the amplitude. The 
boundary condition for the axial velocity that corresponds to the time- 
varying bulk velocity is obtained by scaling measured data for a fully 
developed pipe flow so that the calculated flow field may reach a fully 
developed state at a shorter downstream location from the inlet boundary. 
The numerical results presented herein are obtained at 4 diameters 
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upstream of the exit boundary so that the numerical results are not 
obscured by the exit boundary condition either. Numerical tests show that 
the calculated results become independent of the time-step size for 
AT<T/200, where T is the period. The numerical results presented herein 
are are obtained using AT=T/300. 

The calculated time-averaged wall shearing stresses for 

Re= 15000-60000 and the ensemble-averaged wall shearing stresses for a 
few representative cases are shown in figures 8 and 9, respectively. It can 
be seen in figure 8 that the calculated wall shearing stresses are in very 
good agreement with the measured data for the wide range of Reynolds 
number. The ensemble-averaged wall shearing stress for Re = 15400, 

I u b l^ u b = 0-l» and f = 0.625 Hz is shown in figure 9(a). This pulsating flow 
represents an extreme case of a low Reynolds number flow subjected to a 
high frequency oscillation. In such a case, the spatial variation of the 
oscillating velocity is mostly confined in the near-wall region so that the 
amplitude of the wall shearing stress is a several times larger than that of 
the ensemble-averaged centerline velocity (Mao & Hanrathy 1986). It can 
be seen in the figure that the calculated phase difference is in good 
agreement with the measured data while the calculated amplitude is larger 
than the measured data. The over-predicted amplitude is caused by a 
slightly over-predicted turbulent viscosity in the core region of the pipe. 
The case shown in figure 9(c) represents an opposite extreme case of a high 
Reynolds number flow subjected to a low frequency oscillation. The 
Reynolds number is 60000, |ujj|fa b = 0.05, and f=0.325 Hz. It can be seen in 
the figure that the calculated and the measured amplitudes are by far 
smaller than those shown in figure 9(a). The highly decreased amplitude is 
caused by the large Reynolds number and the small amplitude of the 
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ensemble-averaged centerline velocity. The calculated and the measured 
wall shearing stresses for an intermediate case of Re=15400, |uj(/u 5 = 0.05, 

and f = 0.625 Hz are shown in figure 9(b). In each case, the amplitude of 
the wall shearing stress is over-predicted while the calculated shape and 
phase lead of the ensemble-averaged wall shearing stress are in very good 
agreement with the measured data. It can be found in Kebede et al. (1985) 
that the shape, phase lead and the amplitude of the ensemble-averaged 
wall shearing stresses obtained using a k-e turbulence model and a RSM do 

not compare favorably with the measured data. The significantly improved 
numerical results obtained using the LMS indicate that the turbulence 
equations can correctly describe the influence of the Reynolds number and 
the imposed oscillation frequency on the fluid flow. 

The calculated phase lead of the wall shearing stress with respect to the 
ensemble-averaged centerline velocity is shown in figure (10), where <]>=cot 

is the phase angle. It is shown in the figure that the calculated phase lead 
for co + >0.045 (or co + /15>0.003) obtained using the LMS is in better 
agreement with the measured data than that obtained using the RDT which 
is developed to accurately solve pulsating pipe flows at high dimensionless 
frequency. It can be found in Mao & Hanrathy (1986) that the present 
numerical result is in as good agreement with the measured data as that 
obtained using an optimized relaxation turbulence model. It can also be 
seen in the figure that the LMS correctly predicts the rapidly varying 
phase lead that occurs for co + <0.015 (or co + /15<0.001) while RDT and various 
other turbulence models fail to predict such a phenomenon (Mao & 
Hanrathy 1986; Mankbadi & Liu 1992). 

The calculated ratio of the normalized amplitude of the wall shearing 
stress to that of the centerline velocity is shown in figure 11. It can be seen 
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in the figure that the present numerical result for high dimensionless 
frequency range is in as good agreement with the measured data as that 
obtained using the RDT. It can also be seen in the figure that the LMS 
yields qualitatively correct numerical result for co + <0.015 (or co + /15<0.001) 
while the RDT completely fail to predict the experimentally observed 
behavior of the pulsating pipe flow. 

4. Conclusions and discussion 

A low Reynolds number multiple-time-scale turbulence model that can 
resolve the inequilibrium turbulence phenomena down to the viscous 
sublayer is presented. The capability to resolve the inequilibrium 
turbulence phenomena inside viscous sublayer enables the LMS to yield a 
near-wall dissipation rate that is in correct agreement with the DNS result. 
It is shown that the near-wall turbulence structure for the low Reynolds 
number fully developed channel flows obtained using the LMS is in better 
agreement with the DNS results than those obtained using a Reynolds 
stress model. Each term in the budget of the turbulent kinetic energy 
obtained using the LMS is in excellent agreement with that obtained from 
the DNS, and the excellent agreement is caused by the near-wall dissipation 
rate that is in correct agreement with the DNS result. 

Calculations of unsteady turbulent flows using algebraic turbulence 
models, k-e turbulence models and Reynolds stress turbulence models 
show that the numerical results obtained using the algebraic turbulence 
models are in much better agreement with the measured data than those 
obtained using k-e turbulence models or Reynolds stress turbulence models 
(Kebede et al. 1985; Mankbadi & Liu 1992; Mao & Hanrathy 1986). The 
improved numerical results obtained using the algebraic turbulence models 
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clearly indicate that the algebraic turbulence models can describe the 
transient near-wall turbulent viscosity more accurately than k-e 

turbulence models or Reynolds stress turbulence models. However, 
optimized algebraic turbulence models can not resolve the turbulence field 
of various complex turbulent flows and lack a theoretical background (Mao 
& Hanrathy 1986). It is shown that the low Reynolds number multiple- 
time-scale turbulence model (LMS) yields highly improved numerical 
result for the pulsating pipe flows at a wide rage of the dimensionless 
frequency than the algebraic turbulence model derived from the rapid 
distortion theory (RDT). For pulsating pipe flows at high dimensionless 
frequency, the spatial variation of oscillating velocity is mostly confined in 
the near-wall region. Therefore, a turbulence model that yields only a 
slightly deteriorated numerical results for steady flows may yield a largely 
deteriorated numerical result for pulsating flows. It is shown in the fully 
developed channel flow calculations that the near-wall dissipation rate that 
attains its peak value at the wall yields improved near-wall turbulence 
structure than the near-wall dissipation rate that attains its peak value at 
y + = 12. The improved numerical results for the pulsating pipe flows 
obtained using the LMS is attributed to its capability to resolve the near- 
wall turbulence structure more accurately than other turbulence models. 

Numerical investigations of various complex turbulent flows show that 
the multiple-time-scale turbulence model yields as accurate numerical 
results as those obtained using an optimized turbulence model for each 
flow (Kim & Chen 1989; Kim 1990, 1991; Kim & Benson 1992). These 
accurate numerical results indicate that the developments of the fluid flow 
and the turbulence field in complex turbulent flows depend strongly on the 
inequilibrium turbulence. The improved numerical results for the fully 
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developed channel flows and the pulsating pipe flows obtained using the 
LMS indicate that the development of the near-wall turbulence field also 
depends strongly on the inequilibrium turbulence. The capability of the 
multiple-time-scale turbulence model to resolve various complex 
turbulence fields is attributed to its capability to model the cascade of 
turbulent kinetic energy and to correctly resolve the inequilibrium 
turbulence phenomena. 
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Figure 3.— Turbulent flows in inequillbrium state (Kim & Benson 1992). 

(a) Spectral density curves, k: wave number, E Energy spectral 
density, (P/€ t )A> (P/c t )B > (P/ctJC, 

(b) Semi-emperica! data for eddy viscosity coefficient. 
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Figure 4. — Fully developed channel flow at Re(uJ = 180* 

(a) Mean-velocity profile. 

(b) Turbulent kinetic energy profile. 

(c) Reynolds Stress. 

(d) Dissipation rate. 
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Figure S. — Ensemble-averaged velocity gradient at the wall for 
pulsating pipe flows. 

(a) Re « 1 5400, f = 0.625 Hz, |u b |/u b - 0.1 . 

(b) Re = 1 5400, f = 0.625 Hz, |"u b |/n b = 0.05. 

(c) Re = 60000, f - 0.325 Hz, | u b j/u b » 0.05. 
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Figure 10.— Phase lead of velocity gradient at the wall. 
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Figure 1 1 .—Ratio of amplitude of velocity gradient at the wall to amplitude of 
centerline velocity. 
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